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Chapter 1 
INTRODUCTION 

There are raany transport processes which occur in 
nature and in man-made devices in which flow arises simply 
due to the gradientsof density^ temparaturey^ and/or chemical 
composition in a body force field, such as the gravitational 
field* Ever since the pioneering efforts by Lorenz fl] in 1331/ 
such processes have been of considerable interest to engineers 
and scientists because of their numerous applications* 

Processes in which buoyancy as the driving force arises 
solely due to the temperature difference have received consi- 
derable attention for both steady and transient, internal and 
external, and laminar and turbulent flows with several 
additional conditions and effects such as combined free and 
forced convection, etc* However, buoyancy effects resulting 
from concentration gradients in multicomponent mixtures can 
be just as important in generating fluid motion as the 
t< 3 mperature gradients, as pointed out by Gebhart and Pera [ 2 ]* 
Fields of interest in which combined heat and mass transfer, 
under the condition of free convection are frequently 
encountered are : the evaporation of water from the surface 
of a water body in the absence of strong winds, as from ponds 
and lakes* drying processes in nature; distribution of 
tonparature and moisture over agricultural fields and groves 
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of fruit trees/ damage of crops due to freezing/ formation 
and dispersion of fog/ pollution of the environment/ technological 
applications such as design of chemical piocossing equipment/ 
t^tc* 

In a largo number of important applications# however# the 
convective process is neither predominantly natural nor 
predominantly forced/ both modes being significant# The question 
then is whether the forced convection masks the natural convect- 
ion# and if so# under what conditions? The answer to this 
question is provided here# 

Consideration of transient natural convection is also 
important in many technological applications# since the heat 
transfer rates vary considerably during the transient stage# 
lor given energy inputs this may result in over-heating and 
in consequent damage to various components of the systai:\s, 
furnaces# electronic systems etc# which have# therefore to be 
designed to withstand the transients during the start up and 
shut down operations# We therefore consider the unsteady 
na'tural and forced convection in the presence of tanperature 
and concentration gradients over a vertical flat plate# The 
flat plate provides the simplest geometry so that effects 
other than geometrical can be isolated# 

Earlier developm c^nt s : 

One of the earliest studies with canbined heat and mass 
transfer known to us is that of Scaners [3]# It is concerned 
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with the combined thermal and species diffusion driven flow 
that would arise adjacent to a wetted isothermal vertical 
surface in a non-saturated atmosphere. The condition of very 
small diffusing species concentration was used and an integral 
method analysis vras carried out for uniform surface temperature 
and uniform diffusing species concentration* The principal 
results were a transport relation and the indication that a 
combined driving force might be written in which the species 
diffusion contribution is modified by the square root of the 
Lewis number/ i*e/ Le. The analysis is expected to be 
reasonable around Prandtl and Schmidt numbers of 1*0/ with one 
buoyancy effect being very small compared to the other* 

Mathers et al* [ 4 ] formulated the same problem in terms of 
the boundary layer differential equations resulting from 
momentum/ energy/ and chemical species conservations at 
low concentration* Neglecting inertia effects the resulting 
equations were solved on an analogue computer for Pr = 1.0 
and Sc = 0*5 - 10*0 for ratio of species and thermal diffusion 
buoyancy effects of 1.0 and 0*5* The resultant transport 
information appears to support the VXe factor of Somers. 

Possibility of similarity solutions for combinad buoyancy 
effect flows formulated within the limitations assumed by 
Gill et al. [5] were considered by Lowell and Adams [s] • The 
only similarity solution found was that on an isothermal 
vortical surface. Results of numerical analysis of above 
similarity formulation were presented for a subliming organic 
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surface in air* ISfumerical difficulties that arose when buoyancy 
effects were opposed were noted and termed as flow instabilities/ 
without any satisfactory justification. 

Den Bounter [ 7 ] reports an experimental study of 
simultaneous thermal and chemical species diffusion by an 
elactrochamical mesthod between a vertical copper plate 
maintained at constant temperature and a copper sulphate- 
sulphuric acid solution* Measurements were made with the two 
buoyancy effects aiding and opposing each other. The mass and 
heat transfer parameters/ correlated in terms of a combined 
buoyancy effect/ calculated with the VUe term/ agree well with 
a singlvs curve for the effects aiding each other* However for 
the two buoyancy effects opposing each other the disagreement 
from the single curve is random and over 30 percent in magnitude* 

Bottemanne [s ] has also considered steady state simulta- 
neous heat and mass transfer along a vertical flat plate* Solution 
to the boundary layer equations was obtained only for Pr = 0*7l 
and Sc = 0*6 3. His theoretical solution agrees well with his 
ejg^eriments on heat transfer with simultaneous water evaporation 
into air* 

The problam of combined forced and natural convection 
has been treated by Lloyd and Sparrow [9]/ covering conditions 
ranging from pure forced convection flow to combined flows with 
strong natural convection contribution/ for an isothermal 
vertical flat plats* But buoyancy force arising only due to a 
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temperature difference is considered# The method of similarity 
is employed and numerical results are presented for Pr values 
ranging from 0»003 to 100# 

Some ejqperiraontal worl< has also been carried out on 
mixed convection from a vertical surface tut n;jglecting 
buoyancy effect duo to concentration gradient# Klie-gel [lo] 
employed interferometric methods to determine the local heat 
transfer rates from a vortical surface located in am air 
stream# This data was found to be in very good agreonent with 
the analytical results of Lloyd and Sparrow# 

Gryzagoridis [ll] considered the natural convection flow 
over an isothermal vertical surface X'j-ith aiding external flow 
and determined the local velocity and temperature profiles and 
the heat transfer rates# Good correlation between theory and 
e>qoeriment was obtained for the tsnperature profiles# The 
correlation was found to be very good for the heat transfer 
results/ as ej?pected from the agreement of temperature profiles* 
The limits of forced and free convection regions wore also 
determined for Pr = 0#72# 

Numerical solution for boundary layer equations for a 
transient free convection (buoyancy effects due to temperatiare 
gradivcnt alone) over a vertical surface /subjected to a. step 
change in the surface temperature/ have been obtained by Heliums 
and Churchill [ 1 2 ] • The results converge to steady state values 
at large time and show a minimum in the Nusselt number during 
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the transient stage^ as found earlier by Siegel* 

Gebhart and Pera [ 2 ] studied laminar natural convection 
flows resulting from combined buoyancy mechanisms over a 
vertical flat plate in terms of similarity solutions* Over a 
range of Schmidt numbers, both aiding and opposing buoyancy 
effects were considered for air and water and solutions were 
obtained* The results show many interesting effects on velocity, 
heat and mass transfer, and on laminar stability* In this 
study stratification is neglected and only power law variations, 
t~t^ = N^x^, c - c^ = x^ are considered in order to make the 
similarity solution work* Gebhart and Pera studied the problem 
for steady state* However, the added constraint, brought in 
because of the transient terms, does not allow problems of 
practical Importance to be studied by the similarity variable 
method, its application being restricted to particular forms 
of the surface temperature distribution* A numerical integration 
method may be considered but the procedure is complex one and 
does not yield desired physical insight into the process* 

More recently Callhan and Marner [l3 ] have considered 
transient laminar free convection along a vertical, isothermal 
flat plate, arising from buoyancy forces created by both 
tanperature and concentration gradients* The coupled nonlinear 
partial differential equations are solved numerically using 
an ejcplicit finite difference scheme* Results were obtained for 
Pr = 1*0 and range of Schmidt numbers and for aiding mass 



diffusion buoyancy forces (N > 0)#Steady sta te local Nusselt 
and Sherwood numbers were compared with the results of Gebhart 
and Pera. Largest deviation was observed at X = 0*10 with a 
difference of 4*2 percent while excellent agresnent was found 
at the trailing edge of the plate with a difference of only 
0#75 percent. 

Present Work : 

The purpose of the present study is to investigate the 
problon of transient, laminar, combined forced and natural 
convection along an isothermal vertical plate which is subjected 
to a step —change in tanperature and concentration# The study 
covers conditions ranging firom pure natural convection to strong 
forced convection# The coupled nonlinear partial differential 
equations are solved numerically by a highly implicit finite 
difference procedure# 

Thera are many interesting aspects of such flows, such 
as the resulting transport characteristics, the influence of 
combined buoyancy force effects and combined free and forced 
convection on the stability of the boundary layers, and the effect 
of the values of relative transport parameters, the Prandtl 
and Schmidt numbers# Of particular interest in this study are : 

i) The effect of the buoyancy forces due to' mass transfer 

on the transient velocity profiles, temperature profiles, 
Nusselt number and Sherwood number, 
ii) The effect of free stream velocity on the Nusselt and 
Sherwood numbers, and on the transient velocity, 
temperature and concentration profiles# 



Chapter 2 
ANALYSIS 

-•1 Governing Equations : 

Fundamental physical processes that occur in natural- 
convection flows are essentially the same as those occuring 
in any fluid flow and diffusion processes* Therefore the 
basic equations used to interpret and analyse natural— convection 
flows are the same# There is however one fundamental difference 
between natural convection and forced convection* In natural 
convection fluid motion arises mainly from buoyancy and not 
from imposed motion or pressure difference# The tuoyancy force 
arises due to the action of body force, usually gravity, on 
the density differences in a body of fluid, which result;^ from 
temperature (and/or species-concentration) differences, which 
in turn are governed by the type of diffusion processs present# 
The diffusion processes which may be occuring simultaneously are 
coupled together resulting in much greater complexity and 
difficulty in treating this problsn# 

The basic i^quations are continuity, Navier-Stohes, energy 
and mass diffusion resulting from various conservation laws# 
These equations, j-n general form are [”14] 

continuity ; 


II + V,(PV) = 0 


( 2 . 1 ) 
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momentum s 

= pg - \rp + jii.v^'V + ^ v(v*v) (2»2) 

energy : 

PSf =p[||^ CV.;)e] 

= V>(KVT) + q"'- p V,V + ll(p (2*3) 

species : 

^ + (v.v) c = V.(DVc) + c'" (2.4) 

g V- 

where V is the fluid velocity vector, T is the temperature, e 
is the specific internal energy of the mixture, c is the 
concentration of a single diffusing species defined as the 
ratio of mass of the species (in a given volume) to the mass 
of mixture in the same volume, and c*" are the rates of 

energy and species generation respectively per unit volume, 
is the gravitational force per unit volume, p is the pressure, 
M, K, D are the molecular transport properties namely dynamic 
viscosity, conductivity, and mass 'diffusivity and <P is the 
function associated with the dissipation of energy. 

2*2. App ro xlmat ion s : 

Considering the fluid to behave as a perfect gas, we can 
rewrite the energy equation (2*3) as 

jyp ^ ^ Dp 

pGp ^ = V.(KVT) -h q"' + ^ + 


( 2.5) 
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where 6 is the specific heat at constant pressure for the 

Ir 

mixture* 

The principal difficulties in the above equations ( 2*1) * 
(2»2)j (2*4) and (2*5) result mainly from the possible 
variation of transport properties fi ^ K and D on the one hand 
and density p on the other hand* Since /I / ^ and D are 
dependent primarily on temperature, an appreciable variation 
occurs only in those processes involving large temperature 
differences* Hence these properties are assumed constant here* 
Their variation can however be easily accounted for in the 
numerical method* 

The density differences are approximated/ for processes 
not involving large tamperature differences, by the Boussinesq 
approximation [is]* This simplification renders the continuity 
equation to the constant density form, and introduces into 
the momentum equation a buoyancy force arising from both 
temperature and concentration differences* 

For natural convection flows, frcwn the hydrostatic 
considerations*the pressure gradient vP in the remote ambient 
fluid is Poo g where is the density of ambient fluid* 

p g - VP = g ( p - 1^) • 

How since the surface is vertical and Go~ordinate x is assumed 

positive upwards, as shown in fig* 2*1 the only term of the 

3P 

body force is - p g„ and ~ = -p - g 

* variation with concentration neglected <&ie to low mass fraction 
of diffusing species* 




Body force 
(Gravity field) 


Fig. 2.1 Coordinate system for combined force 
and forced convection flow over a 
vertical flat plate. 
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/. P g - VP = - P) ( 2 . 6 a) 

The series eixpansion of (P^ - P) in terms of T, p and c, at a 
given location can be written as 

(P^ “ P) =P3 (T-TJ + pg'^c-cj (2.6b) 

x-zhere 3 is the volumetric coefficient of thermal expansion, 3 '''^ 
is the volumetric coefficient of expansion with concentration, 

T^ and Cj^ are the temperature and concentration respectively 
in the free stream. 

DP 

The pressure term in the energy equation ( 2*5) is 
negligible for gas flows of small vertical extent as is the 
case here. The viscous dissipation term is also negligible 
for small velocity flows. Moreover the rate of energy 
generation g"' and the rate of species generation c'" due to 

chemical reaction are assumed to be zero. 

2.3* Governing Equations for Natural Convection ; 

With all approximations and assumptions discussed above 
the equations become 


v.v = 0 

( 2 * 7 ) 

P[|~ + (V.v)v] = Pg^(T-Tj + 

p ic-cj + MV^v 


, ( 2.S ) 

P Cp [|| + (v.v)t] = ,k t 

(2.9) 

1^ + (v.v)c = D c 

( 2 . 10 ) 
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The present study is concerned with the simple case 
of two dimensional flow where V = (u,v)» In terms of the 
co-ordinate system shown in fig# ?#1/ the equations ( 2#7) to 
( 2.10) can be written as. 


continuity : 


3u ^ 
3x 3y 


( 2#1 1 ) 


x-momentxim : 


P(|^ + u + V ||) = P g|3(T“Jr^)+ Pg3*(c-c^) 

+ II 

3x 3y 


v^omentum : 


( 2 » 1 2a) 


energy 


jpecies 




PP . 3*^ . — 3Tn 3%* , 9 


fl + - fi + - If = ^ 


( 2.1 2b) 


( 2.1 3) 


( 2-14) 


These equations can be further simplified by carrying out an 
order of magnitude analysis [is]# Let us first non-dim an s ion alize 
the variables as 


- U “ V - X 

u=^,v=g,x = x: 


/ 7 = 2/ 6 


OO CO 

/ C = *2-— — 

\sj OO OO 


where L* U are the characteristic length and velocity, T and 
are the wall temperature and concentration# Equations (2»11) 


to (2#l4) in dimensionless form are 
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du ^ 3v 

3x ay 

od) o(l) 


o 


( 2 * 15 ) 


du 

3T 




NC 


a X 


ay 


^ (5^ + ^) 


0 ( 1 ) o(l)o(l) o(6)o(^) 0 ( 1 ) 0 ( 1 ) 


VGr ax 

oCl) 


atr' 


at 


+ 


+ 


2 — 

a v n 

— o' 


ay 

o(-^) 


V .ax = 

ax ay ^/■Gr ax" ay' 

0 ( 6 ) o(l)o( 6 ) o( 6 )o(l) 0 ( 6 ) o 4 ) 

O 


( 2*16a) 

) 

( 2.16b) 


ae 

at 


ax 


9®. 

ay 


i— (ii + 

Pr'/’Gr ax^ ay^ 


( 2.1 7) 


0 ( 1 ) o(l)o(l) o(6)o(-^) 


ac ^ - ac . 


at 


- ac 

V 


ax 


0 ( 1 ) 




ay Sc ax^ 


o(~) 


( 2*13) 


0 ( 1 ) o(l)o(l) o(6}o(4') 


0(1) oC-ir) 


where N = j^*(c^—Coj,)3/[g(T^-<r^)]measures the relative 

importance of chemical and thermal effects in causing the 

density differences which create the natural convection effect. 

Gr = gg Is the thermal Grashof number, 

• * 3 2 

Gr*= gg L d^~Cjjj,)/i' is the mass Grashof number, Pr = 


K 


D 


is the Prandtl number, and Sc = ^ is the Schmidt number, 
and where the order of magnitude of each term is written 
under it according to the following estimation* 

In the continuity equation ^ must be of order 1 since 

ax ^ 

•• •• 3 ^ 

both u and x are of order 1. Therefore is also of order 1* 

. ay 
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Since y is of order 6 for the hydrodynamic boundary layer, 

V is also of order 6 • Clearly 6 is the hydrodynamic boundary 
layer thichness- In the x-direction momentum equation (2»16a)/ 
we can therefore neglect - 0(1) as compared to = 0(~) . 

ay" 6^ 

Also /Viscous forces must be of the same order as the inertia 
forces and as the buoyancy forces* Therefore 6 = 0(Gr » 

Considering the y-momentum equation ( 2 -16b) it is observed that 
the inertia and viscous terms are of order 6 or less- Hence 
as compared to the x-momentum equation tfie whole of y-momentum 
equation can be neglected- 


For conduction and convection terms in the energy 

equation (2-1 7) to be comparable 6^ = thermal boundary layer 

thictoess = 0 (Pr ^ Gr - Similarly an estimate of the 

concentration boundary layer thickness 6_ is 6 = 0(Sc ’^Gr*"^^^) . 

9 9 c c 

Also the terms -2—?? and are negligible as compared to 


a-x^ 


3x' 


3y2 


and 


3 c 


“g respectively, and can therefore be neglected* 


3y 


This analysis of relative order of magnitude yields the 
following system of boundary layer-tyiDe equations governing 
the distributions of U/V/T and c for transient free convection 
over a vertical flat plate- 
continuity ; 

( 2-19) 


3u 3_v 
3x 9y 


0 


momentum : 


3u , „ 3u . __ 3u 


at ^ 9x ‘ ' 3y 




+ V ^ = y — 5 + iBg(T-Tj-f3%(c-cJ 

3y 


( 2 - 20 ) 
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energy : 


8T 


ax ay ^ 


species : 


ac . , ac ac _ p, a^c 

g-c + ^ a3E + ^ ay ^ ^ 


( 2. 27) 


( 2 . 22 ) 


where we have returned to dimensional variables u,V/T and c 
and where v is the kinematic viscosity (Ji/P) and a is the 
thermal diffusivity (1^” PC ) » 

hr 

2*4 Boundary Conditions And Initial Conditions : 

For an isothermal plate^. the boundary and initial 


conditions 

are 





’^^x^y/O) = 

0/ 

T(x,y,0) = 

• 

c( x/y*0) 

= Cco 

u(x/0,t) = 

0/ 

T(xxO*t) = 

/ 

c( x*0,t) 


u(0,y*t) = 


T(0,y,t) = 

• 

# 

c(0*y^t) 

= 

u(x,«>,t) = 


TCx/^^t) = 

• 

# 

c( X, »,t) 

~ ^oo 


( 2.23) 


vCx^y^O) = 0/ 
v(x^0,t) = 0/ 

where u^ is the free stream velocity in the same direction 
as that induced by free convection i»e» aiding flow. When 
Uoo = O the problan changes to that of pure natural convection* 


2»5* Non-dimenslonalisation : 

By the following choice of dimensionless parameters 
suggested by the above order of magnitude analysis 
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UL 

172 


V Gr 
T-^oo 

e = sr-isr- 


V = 


/ c = 


L Gr' 
V L 


y 


t V Gr 


, 1/2 


/ = 


V Gr 


c-c 


c — c 
w «> 


1/4 


equations (2*19) through (2»22) and initial and boundary 
conditions are expressed in dimensionless form as 

du , dV 


■ 53 ?: + a? 


= o 


( 2.24) 


dU rr dU ^ .. 


^ TT i®. 4 - V 

aT ^ ax ^ 


12 + u ig + V 

aT^ ^ ^ 


- ilf® + e 
5Y ~ 2 ^ ® 

aY^ 

+ NC 



(2.25^ 

39 1 3^9 

3Y Pr 




( 2-26) 

3C 1 3^C 

3? = ^ ^ 




(2.27) 

9(X,Y,0) = 0 

* 

9 

C(X,Y,0) 

= 0 


9(X,0,T) - 1 

♦ 

/ 

C(X,0,T) 

= 1 

i 

1 (2.28) 

e(0,Y,T) = 0 

m 

i 

C(0,Y,t) 

= 0 


e(x,«>,T) = 0 

# 

9 

C(X,<»,T) 

= 0 



V(X,Y,0) = O 

v{x,o,t) = O 


where U, 


u L 

00 


= Re Gr 


.“ 1 / 2 . 


y Gr^/^ 

Equations (2.24) through (2.27) and (2*^) show that dependent 
variables U,V,0 and C are functions of the dimaisionless spatial 
coordinates X and Y, dimensionless time T , and dimensionless 
parameters N^Pr^Sc and Note that for forced flow, T =s (t u^)/(LU 

which is a proper dimensionless time. 
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N measures the relative impoartance of chemical amd 

thermal diffusion in causing the density differences, which 

create the natural convection effect* N is equal to zero 

when there is no species diffusion body force and becomes 

infinite for no thermal diffusion* The momentum equation 

( 2 * 25 ) indicates that N is positive or negative according as 

the mass diffusion forces aid or oppose those of thermal 

“ 1/2 

diffusion* = Re Gr * is the parameter which determines 

the extent of contribution of forced convection to natural 
convection* = 0 is the case of pure natural convection 
while > 0 is aiding flow and < O is the opposing flow* 

Negative values of can cause separation [ 16 ] which can not 
be handled by the boundary layer type equations (2*24) to (2*27)* 

We therefore consider > O only* 

In the case of pure forced convection, the Prandtl number 
(Pr = V/a) relates the relative thicknesses of the ra«Knentum 
and thermal boundary layers 6 and 6^ respectively* Similarly 
the Schmidt number (Sc =V/D) in the pure forced convective 
mass transfer, relates the momentum and concentration boundary 
layer thicknesses for 6 and 5 ^* 

However, in the case of free convection or combined free 
and forced convection, in the presence of mass diffusion 
contribution to the buoyancy force, the relationship amongst 6 , 

64 - and 6 _ is extremely conplex and depends upon Sc, Pr, U and 

w 'V-.' ■ 'OO 

the buoyancy ratio parameter N* 

* For N ■* «», Gr should be replaced by Gr* for non-dimensional izatioi 
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2-6* Heat and Mass Transfer Analysis : 


It. is a common practice to express heat transfer and 
mass transfer characteristics in terms of the flux rate divided 
by the temperature or concentration difference, causing the 
heat and mass transfer respectively. This ratio defines the 
heat and mass-transfer coefficients h and hj^ respectively, x^rith 
the help of which we can define the instantanesous local Nusselt 
and Sherwood numbers as 


Nu 

X 


Sh 


X 


hx _ '5 ^ 

K' - 

^ __Ai< 


( 2 » 29 a) 

( 2»29b) 


where q and m are the heat and mass flux rates respectively-. For 
uniform values of tanperature and concentration differences 
(T^^-T^) and (c^^-c^g), local values, h and hjj , and average values 
h and h^^ a.rve of int-erest. 


These values are ejxpressed in the dimensionless form to 
obtain the instantaneous local Nusselt and SherxArood numbers 
respectively as follox-fs 


Nu^ = 


(P) I X 

° Y=0 


Sh. 


„V4 


= - (|§) 1^^ X 


( 2- 30a) 

( 2* 30b) 


is a constant, let us merge it into Nu and Sh 


Since Gr 

so as to redefine the instantaneous local Nusselt and Sherwood 


numbers as 
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= -X (§§> I 


Y=0 


Y=0 


The instantaneous mean Nusselt number is defined as 




where h = J" h dx = / h dX 


This yields 


Nu„ = ; (- 1^) 


d^' 


dX 


( 2- 31a) 
( 2.31b) 


( 2.32) 


( 2.33a) 


O ” ■ Y=0 

Similarly the instantaneous mean Sherxifood number can be found 
from 




( 2.33b) 



Chapter 3 


FINITE DIFFERENCE FORdULATION 

Solutions of the coupled continuity, momentum, energy 
and species equations (2.24) ~ (2-27) subjected to initial and 
boundary conditions (2»23), were obtained using a numerical 
marching procedure^ Numerical marching procedures are methods 
in which the solution is obtained in a step-by-step manner, 
always moving downstream through the flow field, and forward 
in A number of finite difference forms for the represen- 

tation of the above equations are possible* These are the 
iniplicit. Crank— Nicholson form etc* 

The choice of finite difference representation depends 
on many factors, including the problem itself, and the size and 
speed of computation desired* Explicit difference representa- 
tions are those in which the unknown quantities in the equation 
may be solved for one at a time, as each step in the marching 
direction is taken* While sinple to T-ork with, ejsplicit methods 
are prone to instability and require Impracticably small step 
sizes in order to ensure stability* Implicit representations 
require the solution of a set of simultaneous equations for the 
unknowns as each step is taken in the downstream direction or 
in time* They have no stability constraints and thus allow 
large steps to be taken without any problem* 



22 


However, almost all formulations for transient problems, 
appearing in the literature are ej^liclt* Caliban and Mamer 
I] 13 ] too, have solved the present problem using ej^licit 
r presentation. We use the Implicit formulation nevertheless 
in order to take advantage of its universal stability so long 
as U > 0« This means in particular that there is no restriction 
on the size of stps At or AX. In using the eplicit 
formulation, such small values of At and AX are required 
that computations can become extremely time consuming# Since 
the matrices encountered in the implicit formulation discussed 
here are tridiagonal, the coirputational time required for a 
complete set of calculations at each step is approximately the 
same as that required for the explicit method, but much larger 
values of AT smd AX are pennitted by the implicit method. 

Thus the implicit method seems to be superior to the eplicit 
method. 

3 .1 Finite**Diff erence Ecpiations : 

A non-uniform finite difference grid (Fig. 3*1) is 
imposed on the flow field. Stp sizes AX, AY and AT are 
taken repectively in the X and Y direction and in time# 
Subscripts j+1 and k Indicate the location of the point under 
consideration in x and y directions repectively, while the 
subscript (i+1) indicates the current time. 

The difference fora selected for equations (2*25)/ (2*26) 
and (2# 27 ) is highly implicit in that not only are all 
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Y -derivatives evaluated at j+1 but the coefficients of non- 
linear convective terms are also evaluated at j+1 * This is 
essential for free convection flow since if the usual implicit 
form is chosen it results in U velocity profile decreasing 
linearly from the plate to zero at whatever value of Y is 
chosen as infinity# This result is obviously incorrect# Hence 
the highly implicit scheme is used# 

The finite difference forms chosen for equations (2*24) 
to (2# 2?) are 


u 


i +1 f Ic+l t i-rl 
AY 


- V. 


( 3.1) 
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_1_ 

Pr 


( AY) ^ 


( 3-3) 


C . 

J 


— c 

+ 1 /k/ irfl 3 "^1 i ^ 

AT 

^ "^j+l / k+1 , i+1 

" ^ 


^j+1 J k^ i+1 / k, i+1 

I 4- TJ - - -..- .-. - - 

i -Vl / k^ 1+1 

^j+l^k+l^i+l “ ^j+l/k-l/i+l 
2(aY) 

(ay) ^ 


( 3^4) 


o 

Truncation error is of o(AX) and o(aY)'‘ for mornsn-tuiri/ 
energy and concentration equations^ and of o(AX) and o(aY) for 
the continuity equation. Since the differential formulation 
given above is non“linear, none of the techniques for linear 
algebraic equations may be employed- However, one very simple 
and effective iterative technique is used here* 

To start with equations (3*1) to (3-4) are rewritten 
using superscripts to indicate on which iteration that value 
was obtained/ for example is obtained on the (l)"*"^ 

iteration while U ^3-1 ^ v , i is obtained on the (l+l)^^ iteration- 
In such a linearised form equations (3-1) to (3-4) are rewritten 

S.S 

y(l-f-l) _ „(1+1) _ ,il+l) 

j-rl,k, i-hl j,k,i-l-l j + l,k+l,i+l j+l,k,i+l 

— - + 

^x Ay 


0 


(3-5) 
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U ^ . “ U . 

j •'1- 1 ^ k/ i-f 1 j "^*1 / k/ i 


AT 


+ U 


( 1 ) 

j +1 ^ k^ irhl 


u(i+i) „ u 

j-l-l f k, i +1 j f k^ i +1 


AX 


V 


( 1 ) 

j+l,k^ i +1 


yd+i) _u(i+i) 

3 + 1 , k+l , i +1 j+ 1 , k— 1 , i +1 


2(AY) 


y(l+l) _ y(l+l) „( 1 + 1 ) 

. 3 + 1 , k+l, i +1 j+ 1 , k, i +1 ^ j+l,k-l,i+l 

(AY)^ 


+ 


G 


( 1 ) 

j+l,k, i +1 


+ NC 


( 1 ) 

j-f 1 / k^ irhl 


( 3*6) 


- p 

ji"l/k/i+l j+^/^fi 

AT 


, .rd+l) 

^ ^ j+ 1 , k, i +1 


0^1+1^ - P 

3 + 1 , k, i +1 j,k, i +1 


AX 


+ V 


( 1 + 1 ) 
j+l/k, i +1 


(1+1) fid+1) 

j+ 1 , k+l, i+l“”j+l,k-l,i+l 


2(AY) 


gd+l) 

j j+ 1 , k+l, i +1 

- _ . 
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( 1 + 1 ) 

j+l,k. 


i +1 


4 - pd+i) 

^j+l,k-l,i+l 


(AY)^ 


( 3*7) 


pd+i) ^ ^ 

j+ 1 , k, i +1 j+ 1 , k, i 


At 


+ U 


( 1 + 1 ) 
j-i- 1 , k, i +1 


pd+l) _ p 

i+ 1 , k, i +1 j,k, i +1 


AX 


+ V 


( 1 + 1 ) 

j+ 1 , k, i +1 


pd+i) _ pd+i) 

j+ 1 , k+l, i +1 j+l,k-l,i+l 


2(AY) 


p( 1+1 ) p( 1+1 ) 4- I***! ) 

3 + 1 , k+l, i +1 3 + 1 , k, i +1 3 +l,k-l,i+l 


(AY) 2 


_1_ 

Sc 


( 3 *3) 
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To continuity equation (3*5) and conservation equations (3»6) 

f 

to (3*3) for manentxam, energy and cdiemical species respectively 
are written in more useful form as 


V\ 


(1+1) _ ^(1+1) AYf./(l+l) » 

j+l,k+l^i+l “ ^j+l,k,i+l 5x ^ ^j+l,k+l,i+l ^j,k+l,i+r 


( 3-9) 


yU) 


r_ni_ _ +r -i-+ a 

(ay) ^ '>( ay ) j+l/k-l/i+l 


u 


( 1 ) 


2(AY) 


(ay) 


j+l,k, i+1 

+ ] 

AX 


„(1+1) 

J+l/k^ i+1 


V 


( 1 ) 


^ -1 j+l,k,i+l 


(AY) 


2(AY) 


u u ^ u 

j+l,k,i ^J+l,k^i+l ^j,k,i+l 


AT 


+ 6 


j+1 ^ i+l 

(1) 


AX 


j+l/k, i+1 


[ 


-1 


Pr(AY) 


+ N , 

j+1 f k^ i+l 


„( 1+1 ) 

j+1, k, i+1 

■ J ® -i 


U 


( 3-10) 
( 1 + 1 ) 


2(aY) 


- _ j+1, k, i+1 

.+ [ 7 ^ + —:^ + — ] 


j+ip^.l,iH^ LAT - pj-(^Y)2 


AX 


V 


( 1 + 1 ) 


od+l) . r ■ , jd*X,i+l ( 1 + 1 ) 

w.. ... +L : 77+ ^ : J 5+i,k+i,i+i 


j+1 , k, i+1 I- pj«(^y) ^ 


2(AY) 

(1^1) 

® j +1 , k, i ^ j+1 , k, i+1 ® j , k, i+1 


AT 

^( 1 + 1 ) 


( 3 - 1 1 ) 


AX 


r -1 „ j+1, k, i+1 J ^ 

Sc(AY)^ 2(AY) j+l,k-l,i+l L^-r Sc(AY) 


. r 1 . 2 

+ [ + 7"7-:':x7 
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y(l+l) ^(1+1) 

. r -1 , ^(1-1-1) 

j+l/k^i+l Sc(AY)^ 2(AY) j-!-l,k+l,l+l 

P „(1+1) p 

j +1 , X, i j4-l f X, i+1 J / X, i+1 

= -{- ^ (3.12) 


at: ax 

Equations (3*lo) to (3»l2) writtoi for X = l(l)n constitute sets 

( Itl ) 

of n linear algebraic equations in n unXnowns ^ i+1 * 

®j+l^X i+1 ^j+l^X i+1 ^sspectively. Value of n is chosen 

large enough, so that on several points of the grid close to 
X = n the U velocities are essentially that of free stream. 
Solution procedure and steps followed for the solution of the 
finite difference equations (3»9) to (3*12) to be ejq^lained 
1 at er • 


3 » 2 Heat and Mass Transfer Solution : 

In finite difference form equations (2»3la) and ( 2*3lb) 
for the instantaneous local Nusselt and Sherwood numbers can 
b e wr it t en as 


Nu 

X 


Sh 


X 


X 


X 


30. - — 40 . .j + 6h.,T o 

j+1,0 J+1,1 Jtl#2 


2(AY) 

2^j+l,0 '■ ^^3+1 A ^ S‘tl/2 
2 Cay) 


( 3*l3a) 

( 3»l3b) 


where three-point forvjard differences have been used to 

rs 0 A C 

evaluate (g^) and (|^) at the plate. These differences 
an error of o(ay)^'. 


involve 
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For convenience in es^pressing the instantaneous mean 
Nusselt and Sherwood numbers^ let f and I represent the 
folloiving values at any instant 


6 =-(|y) I 


and t 


=-<§?> I 


Y=0 


( 3tl4) 


Y=0 


Then the instantaneous mean Nusselt and Sherwood numbers can 
be written as (cf. equations (2* 33a) and (2p33b)) 


1 

Nu^ = ; § dx ( 3 . 15 a) 

0 

1 

Sh^ = / I 6X ( 3 . 15 b) 

o 

Using Singsongs rule [l?]/ these can be written as 


nu , = 
rn 


AX 


[§ 0 + 4 ^ 1 + 252+4 § 3+...+ 2 


( 3 »16a) 


= 

m 


AX 




m—2 m* 




(3.16b) 


vrhcre m must be an even number. 

3.3* Computational Steps t 

For the solution of f inite~dif f erence equations in 
Section 3.1 and Section 3*2 following procedure is followed. 

Starting from the specified initial conditions at time 
T = 0, we marched in time first and the velocity^ tanperature 
and concentration fields, and mean Nusselt and Sherwood mimbers 
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Wf;ra obtain . ;d at tiras t + At , for each X location^ sta.rting 
from the leading edge and marching do^'-m stream* At a particular 
location X and time T , the iterative technique is applied as 
'.oqo 1 a in ed below » 


a) 


b) 


c) 


The f j.rst iteration is started by guessing values 

TT ^ t ) ytCo) p, C O ) O ) 

j+l,k,i+l' i-s-1 ^ '^■j^-.l.^k^i+l jb-l,k,i+l 

These giuesses are the values at the preceding step 

upstream (like U . , , , . etc.) 

3 / 


for 


( 1 ) 

Taking 1 = 0 in equation (3i»10)_, ^ are 

calculated solving a set of (n) simultanecpjs equations. 

(solution procedure to be explained later.) 

The continuity equation (3.9) aaain at 1 = O is solved 

( 1 ) 

for the transverse velocity component V.,, , , , 
in the stepwise manner/ working outward from the plate 


surface* 

calculated at 1 = 0 from 

equations (3*11) and (3*12) respectively/ solving 
sets of n simultaneous equations (solution procedure 
to be explained later*) 

•'-) The entire procedure in st^s b) to d) is repeated 

( 1+1 ) 

for 1 = 1/2/3*** and so on, until 

agree to within a desired degree of 

J + I / K/ 1-r.l. 

accuracy S/ taken here as 10 * Similar degree of 

accuracy Is .set for and 

i+1* 
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It might be noted that this iterative procedure is a 
composite of Jacobi and Gauss— Siedel iterative techniques as 
extended to nonlinear equations* 

Once iteration is complete at a particular location X 
and time t ^ the instantaneous local Nusselt and SherW&od 
numbers are calculated solving equations (3*l3a) and (3#l3b) 
respectively* 

Now another step AX downstream is taken and the 
process is repreated at the same t* When the solution has been 
carried downstream as far as desired* the instantaneous mean 
Nusselt and Sherwood numbers are calculated solving the 
equations (3*16a) and (3*16b)* Then another time step A"^ is 
taken* and again starting at the leading edge* the solution 
is marched spacewise downstream* The whole process is repeated 
as many times as necessary to determine the steady state 
solution such that U ^ , . and U . , ^ , agree to within the 

J+X / Ky i J-rlyKyX-ri 

desired degree of accuracy e* 

3*4 Solution Procedure : 

At each iteration* the system of (n) slimiltanaous linear 
equations resulting from momentum equation ( 3*10) can be written 


in matrix form as 
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„(l) 

^(1) 





„( 1) 




^2 

Q 2 



P q 




Q- 


a 


( 1 ) 


^( 1 ) „( 1 ) ^( 1 ) 

'^n-l ^n-l ^n-1 


a<l> | 3 ( 1 > 

n 


v(i) 

j-i-ly. kv' i+1 


(AY) 


2(AY) 

( 1 ) 




4 ^' 


u 


1 

^ 4 . — 4 . 


( ay) ^ 


J + 1 / k/ i+ 1 
AX 


1 


(AY) 


v(l) 

j+1, k/ i+1 
2(AY) 


( 1+1 ) 
j+ 1 / 1 / i +1 

( 1 + 1 ) 

j+1 r 2/ i+l 

( 1 + 1 ) 

j+1 f 3/ i+1 


U 


U 


( 1 + 1 ) 

j+1 /n-ly i+lj 

( 1 + 1 ) 

j+l,n/ i+1 


^P. 


( 1 ) 


u u ^ u 

j +1 , k, i j+3 / k, i'!-l j / k/ i+1 


(p 


( 1 ) 




(P 


( 1 ) 


(P 


( 1 ) 

n-1 


n ' 
( 3*17) 


AT 


AX 


,(i) 

j+1, k, 1+1 


+ N , ... 

J-hl ^ k/ l-rl 

Similarly the systons of (a) simultaneous linear equations for 
energy and species respectively resulting frcro equations (3*11) 
and ( 3 * 12 ) respectively are 
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( 3#1S) 
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T-i t-i "H r-i 



%, w-i 


trf 4" tHI Hh ^ 4" 4* 

r4+fH4"'^+ iHI4"*^4' 
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Solution of system •( 3 * 17 ) to ( 3»19 ) by the inversion 
method is highly expensive even on the present-day computers 

O 

since the number of steps required is of order (n)"^ and the 
storage requirement is of order As the matrix of 

coefficients in all the systems (3*17) to (3»19) is tridiagonal 
a subroutine TRIDIA [I 6 ] is used which wo rhs efficiently* 

3*5* Convergence and Relaxation : 

In some cases it is desirable to either overrelax 
or underrelax the iterative procedure in order to speed up 
or slow down the changes from iteration to iteration, in the 
values of dependent variables* Under relaxation is very 
useful for nonlinear problans* It is used to avoid divergence 
in the iterative solution of strongly nonlinear problems* 

In the present problem after preliminary calculations 
it was observed that the iterative procedure does not converge 
unless some underrelaxation is enployed* The following 
relaxation procedure was used* 


In the direct iterative procedure the quantities 

^ 1 T and appearing in the 

momentumy energy and chanical species equations respectively 

, T y: ttCI+I) ,-,(1+1) 

take on the values of Sj+i,)c,i+l' h+l,]c,l+l 

respectively after each iteration* In the relaxation procedure 
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employed the values of U's^ 0's and C' s are modified as 


<— 4- 'X /TT(ltl) _ tt(1) 

j+l^k^i+l ^ ^ j+l,k,i+l ^U^^j-i-l,k^i+l ^j-:-l,k,i-f-l 


) 


<- 


Itl ) 

j+i,k,i+i + ^e^^j+i,k,i+i 


0 ^^^ ■ ■) 


( 3.20) 


c(i) <_ 

j-f 1 ^ k^ i+l 


J+l,k,i+l ■*■ J-+l,k,i+l 


j+ly k, i-rl"^ 


where \y , Xq and ?vq are the underrolaxation factors. The 
factors found useful in the preliminary calculations were 
Xu = Xg = Xq = 0.3* However as we go downstream they were 
increased to 1.0* They depend also on the value of U^» 
Eventhough optimum values were not used effort was made to 
choose values by trial and error such that computational 
time was minimized. 

3*6. Selection of Step Size : 

No quantitative statanents can be made about the 
choice of step sizes. However smaller step sizes are preferred 
in regions of more rapidly changing velocity/ temperature and 
concentration profiles, in order to reduce the error. Thus 
a fine mesh size is required close to the plate surface and 
close to the leading edge* Effects of singularity at the 
leading edge in the boundary layer problem can be confined 
to a small region close to the leading edge by selecting small 
mesh size AX until the profiles smooth out somewhat. 



37 


Since implicit formulation is used^there is no 
restriction on the selection of step sizes AT and AX, from 
the point of view of stability# However, since at small time 
heat is transferred by conduction only and mass is transferred 
by diffusion only, and Nu^^ and Sh^^^ are inversely proportional 
to 'Tt in that region, small time-step sizes are used initially 
and then time— step size is increased to reduce the computational 
time# Following time— step sizes were used* 


for = 0.0 and 1 #0 


fX = 0#05 
AT = 0.2 

and 

for = 10*0 

AT = 0*05 
AT = 0.2 


(0 < T < 0.4) 

(0.4 < T < steady state t) 

( 0 < TC < 0.2) 

(0.2 < T < steady state t) 


In order to reduce the number of equations to be 
solved, which can effect a considerable saving in computer 
time, and to keep round-off error accumulated in solving large 
nuraber of simultaneous equations to a minimum, it is necessary 
to use variable mesh size. Fine mesh size in regions of 
rapidly varying velocities and relatively coarse mesh size in 
regions of slowly varying velocities is called for. Variable 
mesh sizes employed in Y-direction were as follows : 
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(0.4 < T < steady state t) 
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(0.2 < T < steady state t) 


In order to reduce the number of equations to be 
solved, which can effect a considerable saving in computer 
time, and to keep round-off error accumulated in solving large 
nuraber of simultaneous equations to a minimum, it is necessary 
to use variable mesh size. Fine mesh size in regions of 
rapidly varyi-ng velocities and relatively coarse mesh size in 
regions of slowly varying velocities is called for. Variable 
mesh sizes employed in Y-direction were as follows ; 
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for = 0,0 

AY = 0.05 (0 < Y < 0.5) 

AY = 0»l5 (0.5 < Y < 2*0) 

AY = 0,30 ( 2.0 < Y < 3.0) 

AY = 0,50 (3.0 < Y < 17 , 0 ) 

for Uoo = 1.0 

AY = 0.03 (0 < Y < .30 

AY = 0.10 ( ,30 < Y < 2.0) 

ay = 0.25 <^2.0 < Y< 4 , 50 ) 

AY = 0.50 ( 4.50 < Y < 11 . 50 ) 

for = 10.0 

ay = 0.01 (0 < Y < 0.1) 

AY = 0.04 (0-1 < Y < 0 . 5 ) 

ay = 0.10 ( 0.5 < Y < 1 . 5 ) 

ay = 0.50 ( 1.5 < Y < 6 . 0 ) 

The variable mesh technique provides no difficulty 
for forward or backward first order differences of error o(h)/ 
where h is the mesh size. However, when central differences, 
either first or second order, are used as in the transverse 
direction* difficulties arise at the point of mesh size change. 
To alleviate such a difficulty, consider the mesh size change 
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from a smaller step size A'5^ to a larger step size AY2 at 
the Y-location k = p in fig. 3.1* Application of the central 
difference form for the first or second derivative at k = p 
requires the value at ficticious k = q# This value .is found 
by passing a parabola through the values at k = (p“l).,p and 
(p+1) to yield 

°j+l/q “ ^5,+l + 2 ^j+l,p+l^ 


where Q is the dependent variable at some instant t • The 
derivatives at k = p are then c?>prox±mated as 



3 ^ 


k=p 


2(AY1) 

^j+l/q ~ ^ °j+l/P 
__ 


+ Q 


j+l,p-l 


( 3.22a) 


( 3.22b) 


where Q.', - is oiven by equation (3.2l) and is the ratio 
j+i#q ■ ^ 

of small step size to large step size. 


At the points of mesh size change equations (3.22a) and 
(3.22b) are used in the momentum/ energy, and species equations 
and at these points modified equations are rewritten, as they 
involve central differences in transverse direction. 


The difference representation for the continuity 
equation involves only first order foirward differences in 
the transverse direction* Hence no sucih modification in 
equation (3.9) Is required.- It is ensured that the proper 
mesh size is used in the appropriate region. 
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Modified ec[uaiti.ons at points of mesh size change k = p are 


Momentum Ecruat±on : 

/ 

- , .- + u'--*-'' { ) 




AX 


, » ^R T,Cl+l) ,,(l+l) 

^ ^ 3+l,p+l,i+l “ j+l,p-l,i+l 


2(aY) 


(3-23) 


r V^ „(1+1) . , w(l+l) . o _fR_ 

LTHT .r^-i.-i^'i •*• 2U V^j-M,p,i+1 ■*■ 2 ^j+l,p+l,i+l 


- 1 , 1+1 


„ y(l+l) . U^+l) 1 


<aY)' 




which is written in the simplified form as 


- 1 ) 


4^-1 


[■ 


4 r +1 




2( AY) 


(AY)' 

.( 1 ) 


•3 u 


j+l,p-J/i-f-l 


, ^j+i, P/1+1 2(1-%) V.+I,p^i+1 24 

+ 


AX 


2(AY) 


(ay) 


R 1 ,. ( 1+1 ) 

2 J 


j+l,p,i+l 


^( 1 ) 

, r(^) (-=L- + 

■*■ *-4+§o J j + 1, P+1, 1+1 


R 


2(AY) 
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^j+i,p,i ^j+Lp,i+i ^j,p,i+i ,,v 

+ + 


AT 


,+NC*i^^ 


AX 


(3.24) 

Similarly energy and concentration ecuations are written as 


.(1 + (^), (1 


^R+1^ ^ j + l,P,i+l 


Pr (AY) 


2(AY) 


3 e 


( 1 + 1 ) 

j+l,P“l, i+1 


+ r — 

La-t 


2(1-Sd) , 

R j+l^p,i+l j+l/P,i+l (i+i) 

Pr(AY)^ 2(AY) ax j+l#P#i+l 


2i 


2$: 




-1 


^( 1 + 1 ) 

j+l,p,i+l 


R (AY) 


)-| gd+l) 

2(AY) ^ j+i/P+l/i+1 


P TI ^ 1+1 ) r\ 

j+l,P,l ^j+l,p,i+l j/P/i+l 


AT 


(3.25) 


AX 


[ 


-d . (^)) Cl - C^)) 

R 


^+1 j+l,p,i+l (1+1) 




Sc(AY)^ 


2(AY) 

2(l-4o) 


u 


3 ^j+i^p-i,i+i 

(1+1) 


Sc(AY)^ 

2 


R^ ^j+l,pd-H ^j+l.pdtl ^ ( 1 ^ 1 ) 

2(AY) ^ AX j+l/P/i+l 


24. 


1 + 1 ) 


R { -1 


r K ( 




R (Ay) 


j+i,p,i+i^^ ^(1+3) 

_ ^ j j u. ^ 


2 (ay) 


j+l/P+1 f i+1 


,.( 1 + 1 ) 

-j+l,p, i 3+1, 1+1 ^j^p,i+l 
+ ___________ 


AT 


( 3.26) 


Ax 
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A constant mesh size of AX = 0*0 2 was used with 50 
steps in X-direction^ till the upper edge of the plate (X = 1 .0) 
was reached. The computational procedure took slightly over 3 
minutes of CPU time .for = 0.0^ , 3c, = 0*2 and N = 2.0 which 
reduced further with Sc = 2.0 and N = 0.0* For = 1,0 and 
10.0 CPU time required was around 1 minute. This time could 
have been shortened further if a variable mesh sizes were used 
in the X-direction, but in the interest of simplicity in 
computer programming constant mesh size was used in the X*- 
direction. 



Chapter 4 


RESULTS AND DISCUSSION 
4*1 Limiting Chec ks ! 

In order to assess the accuracy of the numerical 
procedure/ several cases were solved for pure natural convection 
= 0.0) and resulting results v/ere compared with those of 
Callahan and Mam er [l3]» Excellent agreement was obtained 
for steady state velocity, temperature and concentration profiles 
at X = 1.0 for Pr ^ 1.0, Sc = 0.7, 7.0 and N = 0.0, 1.0, 2-0. 
Trcinsient Nusselt and Sherwood members for free convection 
(U^ = O.O) were also compared with the results of Callahan 
and Mamer [is] for N s= 0,0, 2*0 and Pr = 1.0 and for various 
values of Sc. Excellent agreement was found over the entire 
time interval. Based on these comparisons, it is felt that 
the present numerical procedure can predict both transient 
and steady-state results quite accurately. Below present 
results for Pr = 0.7, Sc = 0.2 and 2.0, N = 0.0 and 2*0, and 
U_^ = 0,1 and 10. The Pr and Sc values are representative of 
a large number of gases. 

4 » r. Velocity, Temperature and Concentration Profiles : 

Fig. 4*1 shows the typical development of transient 
dimensionless X-component of velocity U at Pr = 0.7, Sc = 0*2 
and N = 2.0 covering conditions ranging fiom almost pure forced 
convection (U^ = 10 .0), combined flow with almost equally strong 
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na-'cural and forced convection contributions (U^ = 1*0);. and pure 
natural convection (U^ = 0«0)« The profiles presented are those 
at the upper edge of the plate i»a» at X = 1.0* Numerical 
values are listed in tables 4#4 to 4»6* For free convection 
case i*e» at = 0.0, it is observed that the velocity increases 
continuously with time until at t ^ 2-40 it reaches a maximum 
value and than it decreases slightly to the steady— state value 
at T 2«B0 (with degree of accuracy of e =10 • The difference 

between the tenporal maximum in the velocity profile and the 
steady-state value^ however, is quite small and is imperceptible 
if shown on the figure- That is why it is not shown- 

The phenomenon of temporal maximum in the velocity 
profile is somewhat su3:p rising- It has been observed and 
discussed by several investigators for the problem of transient 
free convection on a vertical plate in the absence of mass 
transfer- Siegel [is] based on an approximate integral analysis, 
was apparently the first to predict such a behaviour- Later 
analysis by Gebhart [l9]. Heliums and Churchill [l2 ] and Kleppe 
and Marner [20 ] all confirmed the findings of Siegel- Callahan 
and Mamer [iS] predicted such a phenomenon for the more complex 
problem involving simultaneous effect of heat and mass transfer 
in which N and Sc in addition to Pr, are the controlling 
parameters- The maximum velocity apparently occurs when the 
buoyancy forces in the fluid are largest, and it is clear that 
both the magnitade of the maximum velocity and the time at 
which it occurs are functions of these three parameters- However 
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the phenomenon of temporal maximxim is not observed for 
corabined free and forced convection (U^ = l.O) and almost 
pure forced convection (U = lo.o). In these cases forced 

OO 

convection dominates the natural convection* 

The time required to reach steady state velocity 
decreases as is increased. This is due to qradual masking 
of the natural convection by the forced convection. Figs. 4.1 and 
4.2 show that for the same value of N, Pr and Sc the velocity 
boundary thickness decreases with increasing U^. This is also 
expected and is in line with the general behaviour of laminar 
boundary layers. Numerical values are listed in tables 4*1 to 
4 .6 , 

The effect of parameter N on the steady state velocity 
profile, again at X = 1.0, is shown in figs. 4.3 and 4.4 for 
Pr = 0.7 and Sc = 0.2 and 2.0 respectively for values of 
Uoo = 0.0, 1.0 and 10.0. Clearly, the contribution of mass 
diffusion to the buoyancy force increases the maximum velocity 
significantly for low values of both 0.0 and 1.0, though, 

the increase is comparatively less for = 1.0, for values 
of = 10.0 and more (forced flow regime) this increase is 
almost insignificant. 

A comparison of figs. 4.3 and 4.4 shows that the effect 
of the contribution of mass diffusion to the buoyancy force 
decreases for higher Schmidt numbers. This may be attributed 
to the fact that the rate of mass transfer in the fluid, v/hich 
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in turn influences the buoyancy force, decreases as the 
Schmidt number increases# In the forced convection regime, 
however, increase in Schmidt number does not affect the 
maximum velocity as esqjected. Numerical values for the same 
are listed to tables 4 #1 to 4»12* 

Pigs* 4*5 and 4*6 d^ict the development of transient 
dimensionless temperature and concentration profiles at 
X = 1*0 for Pr = 0*7, Sc = 0*2 and 2*0 respectively and N = 2*0 
for pure natural convection case i*e* =s 0*0* Numerical 
values of tgnperature and concentration distributions are 
listed to tables 4*4 and 4*10* The temperature and concentration 
distributions to fig* 4*5, much like the velocity profile in 
fig* 4*1|^ increase to a maxlinum value at T ^ 1*20 and then 
decrease slightly to the steady state value at t ^ 2*30* 

Similarly in fig* 4*6 for Sc = 2*0 the temperature and concen- 
tration distributions increase to a maximum value at T ^ 1*6 
and then decrease slightly to the steady state value at t ~ 2*30* 
Though the maximum temperature and concentration profiles for 
Sc = 0*2 occur sooner than that for Sc = 2*0, the time at v^ich 
the respective temperature and concentration profiles reach their 
steady state value is the same* Ej^erimental data obtained by 
Goldstein R*J. and Eckert [21] and Kiel for a vertical 

flat plate subjected to a step change to heat flux also verified 
this interesting overshoot phenomenon* As in the case of 
velocity profile the parameters Sc, N and Pr influence the 

* This does not imply that their transient .behaviour is identical* 
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extent of overshoot in the temperature and concentration 
profiles and the Instant at whidh this maximum occurs* However 
this overshoot phenomenon is not observed when forced convection 
in addition, to natural convection comes into play* 

A comparison of the transient concentration profiles in 
figs* 4*5 and 4*6 shows that the concentration boundary layer 
for pure natural convection (U^ = O.O) is considerably thinner 
for Sc = 2*0 than that for Sc ^ 0*2* But for s= 10*0 i*e* 
when forced convection is dominant, this difference in the 
concentration boundary layer thickness at Sc = 0*2 and 2*0 is 
much less as seen from fig* 4*7 and table 4*6* Similar 
conclusions can be drawn for N = 0*0 by comparing figs* 4*9 and 
4*10. 

Figs* 4*8 and 4*9 show the steady state temperature and 
concentration profiles at X = 1*0 for Pr = 0*7# Sc = 2*0 and 
N = 2*0 and 0*0 re^ectively for all three values of FxxDm 

figs* 4*5/ 4*6, 4*8 and 4*9 comparing transient and steady state 
temperature and concentration profiles for various values of 
N, Sc and it is observed that the concentration boundary 
layer is thicker than thermal boundary layer for Sc = 0*2 while 
opposite is true for Sc = 2*0* This difference in thiclaiess 
decreases as increases* Fig. 4*10 shows the steady state 
concentration profiles at X ss 1*0, Pr »* 0*7 and Sc = 0*2 with 
N and as parameters* Note that the effect of N on the 
concentration profile decreases as is increased* At = 10*0, 

* The correct parameter is the Lewis number, Le = Scy'Pr, tut Pr is 
held constant here* 
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both values of N yield almost the same concentration profile* 
This is e3(p acted since Ixioyancy forces are insignificant at 
such a high forced flow condition* Numerical values are listed 
in tables4*l to 4 *12* 

4*3* Nusselt and Shejrwood Numbers : 

Transient mean Nusselt and Sherwood numbers are shoxvn 
in figs. 4.11^ 4*1 2 and 4.13 for Pr = 0*7^ Sc = 0*2 and 2*0, 

N = 0*0 and 2*0, and = 0.0, 1.0 and 10.0* Numerical values 
are listed in tables 4.13 and 4*14* Initially for any set of 
parameters values of mean Nusselt number and Sherwood number 
are high but they drop drastically with time and approach 
steady state values. During this initial transient has 
negligible influence on both Nu^ and Sh^^ due to the fact that 
heat is transferred by conduction only and mass is transferred 
by diffusion only during this regime. As the buoyancy forces 
due to mass transfer and thermal convection increase, the 
velocity increases sufficiently for N and to influence the 
solution. Because of the overshoot phenomenon observed in the 
t exaperature and concentration profiles (i.e. these profiles 
reaching a maximum before steady-state conditions are reached) , 
a transient minimum is observed in both the Nusselt and Sherwood 
numbers, for pure natural convection. The difference between 
the temporal minimxam and the steady state value, however, is 
quite small, and in most cases is nearly Imperceptible on the 
figures. (See the inset in fig. 4.12). However for = 1.0 and 
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10*0 where forced convection is predominant such a teniporal 
m3.niimim is not observed since the overshoot phencmenon is 
not present* 

Both the mean Nusselt and the mean Sherwood numbers 
show a slight dependence on Sc and as far as the time 
required to reach steady state conditions is concerned, i.e*, 
as Sc increases and as decreases the time decreases* Kowever, 
the time required to reach steady— state conditions is virtually 
insensitive to the parameter N* 

As the Schmidt number increases from 0*2 to 2«0 with 
Prandtl no* fixed at 0*7 and N fixed at 2*0, Sherwood number 
increases substantially but the Nusselt number is only moderately 
affected* 

, Prom figs* 4*14 and 4*15 it is noted that an increase 
in N results in higher Nusselt and Sherwood numbers* That is, 
an increase in the buoyancy force due to mass transfer results 
in an increase in the rates of heat and mass transfer* The 
effect of parameter N on the Nusselt number becomes less 
pronounced with increasing Schmidt number hut its effect 
on Sherwood number is opposite and to a small extent* This 
behaviour can be ejgjlained as follows* It was observed that 
an increase in Schmidt number decreases the concentration 
boundary layer thickness* However, the thermal boundary 
layer thickness is relatively less sensitive to an 
increase in the Schmidt number for fixed values of Pr and 
aventhough the conservation equations are coupled* Thus with 
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an increase in Schmidt number the concentration boundary 
layer thickness becomes thinner than the thermal boundary 
layer thickness# Hence the influence of parameter N on I'TUj^ 
reduces with increasing, Schmidt niimbar, as the buoyancy effect 
due to mass transfer are diminished in the thermal boundary 
layer# On the contrary thermal buoyancy effects become less 
important as compared to buoyancy effects due to mass transfer 
in the thinner concentration boundary layer. Hence the 

influence of N on Sh increases with increased Schmidt number. 

m 
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Fig. 4.1 Transient velocity profiles at X = 1.0 for Pr = 0.7, Sc=0.2,N = 2.0 
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0.02 31 

O.Q03S 

0.9937 
v) , 9 d 9 6 
0 .3891 
o . 3 9 04 
0,0922 
0.09 40 
0,0954 


6 . 201 4 

3.2596 
1 . 47*^0 
2.1130 
1 ,9873 
1 .7377 
1,6303 
1 .'^493 
1.2296 

1 ,'. >8 6 3 
1 .00 32 
3. '">5 2 7 
:).9213 
o .50^8 
0.9 9 DO 
0,0829 


II =0.3 


.60 
.BC 
1.00 
1.20 
1,40 
1 . D 6 

i.Bo 

2.00 
2.20 
2.4 9 

2.00 
2 . 3 0 


%zid 


l.uo 
1.20 
1.40 
1 .00 
1,60 


m^ss6hs: 

' . ri 

3.5796 
1 .9344 
1.4701 
i.?ol3 
1.1297 
1 , <■ * ^ i 4 
0.«779 
n.Q305 
0 . 7 5 F 6 
n.6d01 
0.6 542 
0.6 37 9 
0.6317 

n . 6 iC 4 

o, '=^.3-:9 
0.6818 
0.6327 
2,6332 
"i , * 3 6 5 
0.6337 


3.7191 
1 .Oi‘’7 

I . d -i 6 D 

1 . 2 f F 5 

1 . IG^.J 
1.083d 
1 . 0 27 3 
0 , c 3 6 3 
C • 8 1 2 

0.7576 
').72Q7 
0.7175 
0.7125 
0.7 1C ; 
•0.71 f 2 
'■.71 1’ 2 
0 . 1 < ’ 2 
. 7 1 '■> 3 


F.lFKw'Jo, 

•■I ’J . 

0 , u b 7 ^ 
3.3033 
2.5117 
2,1380 
1.91 iQ 
i.b733 
1.0550 
L.578F 
1.290^ 
1,1593 

1 ,U‘9 75 

i.c67« 
1.0662 
1.0536 
1,0549 
l.b570 
1.O590 
1 , 6 L' 4 

I , U 6 I r 


O . 2 * ' 9 

3.2544 
^,•+95 0 

1.9'091 

1.7'^25 


.2695 
. 1 •! O 4 

* ■■ I •♦ ? 
*b72- 
C.o 85 
I.l dF 
ImUf 
. I- ’ 7 1 

1 riO 


3,7913 
2 , 0 6 2 3 
1 , o 93 Q 
i;5444 
1.1701 
I .u70l 
i:U353 
1.U240 
1 ,6206 
1.0195 
1.0193 
1.0192 
1. 3132 


6 , 1 7 F 8 
3.1373 
2.4279 
2.1191 
1 . 5 'J 4 2 
1.3034 
1 .232b 
1.2097 
1.2011 
1.1988 
1,1982 
1.1 931 
1.1980 


=10. 
JS 


3.793/ 
7 . 0 '65i- 

1 . - J 8 e 

1 . 

t .1 o9. 
1.0.C7 7 


. o 0 
i.0'1 
1.20 
1.4} 

1.3 ' 
1 .b>^ 


o . 1 = V.' 4 

3. 1-1 '-5 
2.^934 
4 , 1 2 c 4 
l.b207 

l,37o9 

1 . 2 6 . 

1.239' 

1.23/5 

1.230*- 

1.2347 

1 , 2 3 J 1 


4 = ' I .'} 


TdOle 4,14 ; 


Transient --e-aa 
?r = 0.7, Sc = 2, 


lusseic dOJ 


.. = 2 . U 

nu^cf^rs for 



Chapter 5 


CONCLUSIONS 

A boundary-layer analysis for transient, laminar, 
combined forced and natural convection along an isothemal 
vertical flat plate subjected to a step change in tanperature 
and concentration has been solved by a highly implicit finite- 
difference method* In order to obtain a solution to this 
problem, the coupled governing conservation equations must be 
solved simultaneously. The canputer code in the Appendix 
is based on a non-uniform mesh in the direction normal to the 
plate ( Y-direction) • The parameters of the problem are : 

( i) the buoyancy ratio parameter, N, ( ii) the Prandtl number, Pr, 
( iii) the Schmidt number. Sc, and (iv) the Forced-free convection 
parameter, = Re/Gr^*^^* Results found for various values of 
the above parameters show the following : 

a) During the initial transient period the heat transfer 
is by conduction only, and the mass transfer is by 
diffusion only, even for strcng forced flow* 

b) After the initial conduction— diffusion rcjgiiaa, 
combined buoyancy forces along i-^ith free straan 
velocity (for combined forced and free flow) generate 
the motion. 

c) The transient velocity, temperature and concentration 
profiles for free convection show a temporal maximum 



90 


over their respective steady state values* 
However, this phenomenon of temporal maximum 
is not observed for combined free and forced 
convection* 

d) The time required to reach the steady state 
decreases with increase in forced-free convection 
parameter 

e) Both Nusselt and Sherwood numbers pass through 
a temporal minimum before reaching their steady 
state values* However, it is observed only for 
free convective flow and not for combined free 
and forced convection. 

f) For mass diffusion aiding the flow both mean 
Nusselt and Sherwood numbers are higher than 
those for N = 0.0 (Pure thermal convection) » 

g) Mean Nusselt and Sherwood numbers are higher 
for higher values of 
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C 
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c 
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M 
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o Li i..'* vj a* T'l O' s r a* L- i. o* «■** a i lie "i 
free f o r ced c a n v a c: t i o n p a r a m e t e t' 

I' e 1 B >i t i o i”i t a c t o r 

s t e p s ;i. z e s n Y d :i e c: t i o n 
t T"' a ii i c.'5 T s Ilf <3 .1, i. e r* s t e p i. z e 't o 
1 a r ^ e r s't e p :i, z e i n Y d ;i. r e r: t i o n 
S step size in X direction 
t t i 111 e-*“ steps 
in X direction 
in Y direction 
at chsnsfe of mesh size 


..Y3i^DELY4 


f"! u m D e r 
n u ill b e r* 


of 

of 


s ‘c e p s 
s t e P e 

3ints 


IP?IGbIR 

: >1{ ;Fc t ;{•; t :H H< 1 1 1 1 1 >!< * 1 1 W- 1 M t % 5{; 1 1 ;{? t $ :Ft t ;F; * f 1 1 1 :K t :}: 1 5K t :=!< $ * * t * >f: sK :K * :* * t- 1 * ■$ 


REAL 

NIJX^' 

NUM 




DIMEN 

SION 

A C 39 

) S*' 

Ev 

( 39 ) 

DIHEN 

SION 

ATHET C 


9 ) f P 

DIHEN 

SION 

CON C 

C' 

) 

»ACO 

DIHEN 

SION 

111(3 


5 

1 ) ? T 

■U .!. M E N 

SION 

NUX ( 

51 

} 

pSHK 

OPEN 

(UNI 

T==:2:i. f 

DE 

y 

ICE” 

OPEN 

< UNI 

*r 2 5 i? 

DE 

y 

ICE” 

OPEN 

<UNI 

T*-23j^ 

DE 

y 

Ii::E=“ 


C Re; 


fO in 


' D SK ? F” I i... E = FF F F * C DR ' ) 
'• DSK '■ f F I LE== ' FFCC 1 6 , BAT 
'DSK'’ !'FILE=='NS16.DAT'' ) 

B r'l d p v i n t o u t t he pare m e t e r s 
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F?E AD < 2 :l. 1- $ ) N ? I P p I Q ? I R p PR p p SMALL p SC p R^^ •' i2 
Rl“ AD ( 2 :l. p >iJ ) DELY :l. p DELY2 p DELY3 p DEL Y4 p DT ALi 1 p ^ ■ 


READ(21?!f<)DELXpMpALPpU:ENF _ ^ , Y4 p 

14 R :i: T E ( 2 5 P 1 :l. :l. ) N P M ? D E L X p D E L Y 1 ? D 1::. L Y 2 p D L L Y p i- 

1 p SMALL •• „ T^p 

:i. :L 1 FORMAT ( 2X P ' NO OF STEPS I N Y D I REC 1 1 UN = / - 

:Lp2Xp''N0 of STEPS IN X DIRECTION == 'plip/ 


2?2Xp'STEP SIZE IN X DIRECTION 


' pF4.2p.'' .y,F4.2' 


o,;LXpF4.2p 'S' 


y a\. y\ y \.} i u.. i w .i. * > ' -• i X ^ .-"i Y 

3 p 2X p ' STEP S I ZES I N Y D I RECT I ON ~ ' pi- 4 . p - ’ 


4 p F 4 . 2 p / p 2 X p ' T I M E S T E P 


' pF4.2p2XpF4.2p « p-7,5.J 


than 

r '*> )C 

"" r< ^ 4 O 

/ y F* 4 ^ 


5p 'CONUERGENCE CRITERIA IS EPS I 
kiRITE(2Sp94) PRpSCpRNpUINT- _ _ ^ 

9 4 F 0 R M A T ( 2 X p ' F o r ' ? 3 X p ' P ° i"' 0 1 1 N 0 » ' p !- ^ ^ ' 

1 p ' S c h m i d t N 0 ™ ' p F 4 « 2 p / p 2 X p ' P s v- s m © t & r ?'• 

2 p 2X p ' %. ' p 2X p ' R©/ < Q r ) ttO » S “ ' p 1- L ♦ 2 ) 

14RITE(20p95) 

95 F0RMAT< 1Xp50< IH-) p/) qHFRWOOD NO»') 

WRITE ( 23 p 888) . , y 

888 FORMAT < 4X p ' TAU ' p 5X ? ' MEAN NUbbEL 1 NO » • ? ■ 

WRITEC23p887) 

007 eormA 1 ( 2X p 50 c I H— ) / .j ... , I'i sLed res^'©3L0dl 

Ms lues of 311 the constants reauireo to 

calculated here 

D Y2n=0 . 5./DELY2 p DY 1 ==0 . 5./DEL Y 1 5 DTAU==DT AU 1 Y'XfUiELX 

DYS-^O . 5/nELY3 p IiY4^=0 . 5/Dl::.L 14 ,, j rnv,:::^. 

DELX=1 » /DELX p D:l. :^a:iELY 1 Ji^DLLX p DL^^DLL t .ciUL - 


., o./DTAU 


_.pp.j3YEl-BYM-i./PR 
.'■i s p DYE4='-=DYM4./PK 


D?«?S!o*S?2«,Y2n.VHl=4.0*DYl*OYi;DY»»g-;5/P-<'”'' 

DYM3-~4 * 0^:DY3*DY3 p DYM4=4 « Otli 1 45iiD 1 4 p D 

DYC2:=DYM2.Y'SCpDYCi==DYMl/SC _ .^4 

nYC3==DYM3p'SC.pDYC4==DYM4.Y bC ,t f APHItl ^ ) 

APH I ==DELY 1 /DEL Y2 p EPH I =-U.iELY2/DbLY -.5 ^ SR Hi j- 1 * •' ' 

A F^’ F I I ;i. ~ 2 4 1 A P H 1 1 A P H I / < 1 " 0 + P H 1 > p A P H 1 2 " <•. h ^ ^ p ^ ^ ^ ^ .j ^ y 

A P H 1 3 == 2 <• * ( :i. 4 ~ A !■•' 1 -I I > ri P H 1 ^ 

BPHIl^=2«t;EPHI>KBPHI/< 1 ^OTBPHI ) J BPHI2- 1 /OOPHITl . ) 
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DO 180 K= 

=2sIP 


Y ( K ) “Y ( K- 

■DtDELYl 

80 

CONTINUE 
DO 190 K= 

= IP2sIQ 


Y ( K ) ™Y ( K- 

■1 >-fDELY2 

90 

CONTINUE 
DO 191 K=: 

= IQ2sIF? 

91 

Y ( K ) “Y ( K" 
CONTINUE 

■l)iDELY3 


DO 192 K= 

^IR2yN 


Y ( K ) ==Y ( K - 

l)tDELY4 

92 

CONTINUE 



$ 

t initial and boundary conditions are Generated and initial sdjessss 
for Li's? U'sp Theta's and C's at Is an iteration counter sin the 
$ i t e r a t i v e n r o c e d u r e ( e f •» s e c t i o n 3 * 3 ) a r e d i v e n 
5 K 

DO 10 K==l sN 

AU < K ) ==U I NF s A V ( K ) ~0 . 0 s ATHET < K > =EXP ( - 1 <• 9*Y < K ) ) 

AGON C K ) ==EXP ( - 1 » 9?!i Y ( K ) ) 

P U < K s 1 ) -= Li I N F s P I H E T ( K s 1 ) = 0 , 0 s P C 0 N < K s 1 ) := 0 s 0 
DO 10 L“1 sH 

LI 1 ( K s L ) ™U INFs THET 1 < K s L ) « 0 s CON I ( K s L ) ==0 « 0 

10 CONTINUE 
C 

C marchind in time 
C 

125 TAU^^TAUf-DTAU 
WR];TE<5s>}t)TAU 
C 

C for hidher timers relaKation factor - ALP is increased 
C 

IFCTAU.GEsl » )ALP=1.0 

C ms re hind in X - direction (downstream) 

J:rO i5X==0.0 
DO 140 L^=d. sM 
X=O<ils0/DELX 


400 ITR"=ITR'fl 

0 1 0 1 i 1 0 n 't s 0 f T 1 *^! 1 I A G Q H A L. c o e f 'f i c i e n *t m a ‘t r :i. s ii d r i 1*1 't n d ‘Hi i d e 

V0c‘‘‘*bor ”f*or ri*iui!!0*nlvuiTi ac-KU-ai/lor’i (cf ^ 00rf<’ i / /*af‘0 aula Ivt^ci 
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DO 20 K^l fXP 
PR0D^=:AM(K)5f«DY:l, 

A<K)::::--DYM:|,--PR0D 
B ( K ) ™ A U < K ) ;{< .0 E L X i 2 . 0 >K D Y H ;l. -f D T 
C(K)---=---DYM;LfPF?OD 
20 CONTINUE 

A ( I P ) :::: A ( I P ) T A P H I 2 1; C ( I P ) 

B ( I P ) ™ B < I P ) A P H 1 3 1 C ( I P ) 

C<IP)::::C(IP)5i«APHIl 

t 

t step size is chsriiged from smsller step size DELYl to Isrder step 
« size DELY2 
t 

DO 30 K==IP2 5'IQ 
PR0ri==-AyCK)^DY2 
A < K ) =^="-DYM2-“PR0D 
B < K ) ™AU C K ) ;{cDELX+2 « 0>1{DYM2+DT 
C':K)=--DYH2fPR0D 
30 CONTINUE 

A < I Q ) = A < I Q ) f BPH 1 2)fcC (IQ) 

B < I Q ) ==B ( I Q ) f BPH 1 3*C ( I Q ) 

C<IQ)=C(IQ)}i<BPHIl 


step size is chsnded from smsller step size DELY2 to larder step 
size DELY3 


DO 21 K==IQ2fIR 
PR0D==AV(K)tDY3 
A(K)~""DYH3“-PR0D 
B ( K ))=AU < K ) $DELX f 2 * 0.fDYH3iDT 
C<K)--""DYH3-f-PR0D 
21 CONTINUE 

A ( I R ) == A ( I R ) i C P H 1 2 C < I R ) 

B ( I R ) ~B ( I R ) -f CPH 1 3-45C < I R ) 

C<IR)==C(IR):{:CPHI1 

t step size is c handed from smsller step size DELY3 to Isrser step 
size DELY4 


DO 31 K=IR2?N 
PROD-AV ( K ) $DY4 
A ( K ) ^=>-DYH4“PR0D 
B ( !•< ) == A U ( K ) D E 1... X i 2 . 0 % D Y H 4 1 D T 
C ( K ) -=--DYM4'f PROD 
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DO 90 K=-IP2'-IQ 
PR0.l1=y(K)>l?DY2 
A < K ) ™--DY£2-PFv'0ri 
B ( K ) ==U ( K ) J^DIELXf 2 , 0>|criYE2fDT 
C ( K ) ~-DYE2fPR0D 
90 CONTINUE 

A ( I Q ) :==A (IQ) iBPH 1 25KC ( I Q ) 

B ( I Q ) (IQ) f BPH 1 3;lcC ( I Q ) 

c<iq):^h:::(iq)5Kbphii 
t 

$ step size is chansSed from smaller step size DELY2 to larsier step 
size DELY3 

DO 8:l. K==:IQ2i'IF? 

PF?0D“M ( K ) :^cDY3 
A < K ) =>-dye3-pf;;od 
B ( K ) -••■U ( K ) >!<DELX+2 . OfcD YE3+DT 
C(K)----=-“DYE3iPR0D 
81 CONTINUE 

A ( I R ) ==== A ( I R ) ■}• C P H 1 2 )fc C ( I R ) 

B ( I R ) = B < I R ) i C P H 1 3 >K C ( I R ) 

C(IR)™C(IR)*CPHI1 

step size is cnansed 'f rom smaller step sizos DEL j -jj to larser step 
size DELY4 

DO 91 K==IR2pN 
PR0D==U ( K ) ^DY4 
A ( K > ="-DYE4--PR0D 
B ( K > =U ( K ) ;f DELXf 2 . 05F;D YE4f DT 
C<i<)=-DYE4fPR0D 
91 CONTINUE 

DO 100 K=i?N 

D ( K ) ( K ) :i?PTHET ( K ? L ) *DELX-f TFIETI ( K ? L ) :$DT 

100 CONTINUE 

D ( 1 ) < ;i. ) - A ( 1 ) :l:THE0 

C A L. L T R I D I ( A ? B C T H E T ? D ? N ) 

Y 

caiculst ion of Theta's (temperatures) at various ¥-■ locations 
>K is comp' 1 et e 

Y 0 ? 1 e merits of T rtlDJ. AuUNaL coGJf 'f'ic ient mat rr >; arid r isn't fiand siioe 
% y e c t o T f o r s p- & c i & s e & u s t i o n ( c f » e ca n * 3 » 1 9 ) a i" e c a 1 c i..i I a t s d 
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DQ :L50 K-ImI'P 
FR0):i==V(K)S:Dy.1. 

A(K)^=^-riYC;i.-PROD 
B < K ) -U ( K ) :=KDELX f 2 <► OsUDYCl-fBT 
C(K)™-"DYCl-iPROD 
:!.50 CONTINUE 

A( IP):=AC IP)-fAPHI2>icC( IP) 

B ( I P ) B ( I P ) -f A P i“l 1 3 5K C ( I P ) 

C(IP):=CCIP)>lcAPHI:L 

5K 

>K step size is chsn^ed from smeller step size DELYl to larSer step 
% size DELY2 

DQ 160 K=IP2?N 
PR0D=M<K)jUDY2 
A ( K ) ~-DYC2-PR0D 
B < K ) ==U ( K ) >f:DELXi2 . 05{cDYC2+DT 
C ( K ) ==-DYC2t PROD 
160 CONTINUE 

A < I Q ) == A ( I a ) i B P l -l 1 2 * C C I Q ) 

B (IQ) ==B < I Q ) 4 BPH 1 3$C (IQ) 

C( IQ)=:C(IQ)5i<BPHIl 

step size is chsnSed from smeller step size DELY2 to larder step 
size BELY 3 

DO 151 K=IQ2?IR 
PR0D==U(K)^DY3 
A(K)™-"DYC3-PR0D 
B ( K ) ==:U ( K ) :icDELX f 2 . 05f:D YC3fDT 
C<K)==--DYC34PR0D 
151 CONTINUE 

A ( I R ) := A ( I R ) i C P H 1 2 C ( I R ) 

B ( IR ) ( IR ) t CPHI3)!£C ( IR ) 

C( IR)=::C( IR):tCPHIl 

ii step size is ch ended from smaller step size BELTS to larder step 
4; size DELY4 

.-•jv 

DO 161 K:=^IR2?N 
PR0D==^U(K):tDY4 
A ( K ) ■•^^-■DYC4-PR0D 
E ( K ) ( K ) f-BEL-X-f 2 . OfDYC^iDT 
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C<K)::=-DYC4iPR0D 
161 CONTINUE 

DO 170 K==1jN 

D ( K ) ( K ) :>i5PC0N ( K ? L ) >KDELX-f CON I ( K » L > JKDT 

170 CONTINUE 

D(l)==D(l)-A<l);icCONO 
C A L L T R I D I ( A j B ? C s- C 0 N p D i- N ) 

){£ 

Y c 3 1 c u 1 B 't i o n o T C s ( c o n c b n t r 3 1 i on 3 1- v s r i o i..! s T “* 1 o c s t i o n s i s c o ift p 1 s t- 
t whether solutions for U''s?V's!'Thets'syC'3 durinsJ iterative nrocedur 
if: are converj-Sed or not is checked* if wes we move to neKt X- location 
4 Qownstream* if not the V3lues calculated in this iteretion sre the 
'1< sluess values for the nent iteration* under re l3K3t ion is ernploaed 
$ r e f e t ' t o s e c t i o n 3*3 fo r d e 1 3 i 1 s 

CALL MAX ( AU ? Li i- N i- EPS ) 

IF ( EPS *LT, SMALL) GO TO 500 
800 DO 120 K-li-N 

AU < K ) :==AU ( K ) -f ALP^ < U ( K ) -AU ( K ) ) 

AO < K ) ==AV < K ) f ALP:?? ( V ( K ) -AU ( K ) ) 

ATHET C K ) -ATHET < K ) -f ALP>i? ( THET < K ) -ATHET < K ) ) 

ACOM ( K ) -ACON ( K ) +ALP:i? ( CON ( K ) -ACON ( K ) > 

120 CONTINUE 
GO TO 400 

500 CALL MAX ( AU s- 0 f N ? EPS ) 

IF ( EPS *GT* SMALL) GO TO 800 
CALL MAX ( ATHET i- THET i-N HIPS) 

IF < EPS *GT* SMALL) GO TO 800 
CALL MAX < ACON y CON i-Ni- EPS) 

IF ( EPS *GT* SMALL) GO TO 800 

the instantanecsus local Nusselt and Sherwood numbers are calculate 
s o 1 V i n S t h e e o. u a t i o n s (3*133) e n d (3*13 b ) r e s f e c t i v e 1 s 

N U X ( L + 1 ) - ( 3 * 0 - 4 * 0 * T H E T ( 1 ) t T H E T ( 2 ) ) D Y 1 
SHX ( L + 1 ) == ( 3 * 0--4 * 0*C0N ( 1 ) i CON ( 2 ) ) :{:DY1 

t writing out the values of U?V?Thete and C et every Y-- location but 
t at X - 1*0 i* e* at the UFF-er edde of the F-late 

IF (J*NE*50) GO TO 786 
W R I T E ( 2 5 j 9 7 ) X ■> T A U ? A L P ? I T R 

97 F0RMAT(2X? ?F3*1?1X? 'TIME = " ? F4 * 2 ? IX ? ■' RELAX FACTOR == 'yFS*! 
IflXy'ITR :::: " ? I 2 !- / ) 
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C ( K ) :--DYC4f PROD 

CONTINUE 

DO ;l.70 K==::i. i-N 

D < K ) ( K ) .KPCON ( K ? L ) $DEI...X-f CONI ( K p L ) 5!?DT 

CONTINUE 

D < :l. ) ^==D < 1 ) --A ( 1 ) 5KCONO 
C A L i... T R I D I C A p B p C p C 0 N p D p N ) 

I c u 1 a t ;i. o n o f C ' s ( c o n c e n t r a t i o ri a t v s r i o u a Y 1 o c a t i o n s i a c o rfi p 1 e t e 
the? solut-i ons for U' spU'spThe?bs''spC‘'s dur-ini=^ i ters’tive procedure 
converged or not is checked « if yes we move to next X-’-locst ion 
nstresiTi« if not the values calculated in this iteration are the 
ss values for the next iteration « under relaxation is siTiPloyed 
er to section 3i-3 for details 


CALL HAXCAUpUpNpEPS) 

IF' (EPS.LTp SHALL) GO TO 500 
DO :l.20 K~:LpN 

AU ( K ) “mAU ( K ) i ALP* ( LI ( K ) -AU < K ) ) 
m ( K ) ™AV < K ) i ALP* ( U ( K ) ••■■Ay C K ) ) 

ATHET ( K ) “ATHET ( K ) -f ALP* < THET ( K ) -ATHET ( K ) ) 

ACON ( K ) ::=AC0N ( K ) -f ALP* ( CON C K ) -ACON ( K ) ) 

CONTINUE 
GO TO 400 

CALL HAX (AUpMpNpEPS) 

IF C EPS «GT. SHALL) GO TO 800 
CALL MAX < ATHET p THET pNp EPS) 

IF ( EPS. GT* SMALL) GO TO 800 
CALL MAX < ACON p CON p N p EPS ) 

1 F ( EPS . GT . SMALL ) 60 TO 800 

e instantaneous local Nusselt and Sherwood numbers ere calculated 
Ivind 'the eeuations (3.13s> and (3.13b) re-spect ively 

NUX ( L-Fl ) ™ ( 3 . 0--4 . 0*THET ( 1 ) -fTHET ( 2 ) ) *DY 1 
SHX ( L-l- 1 ) ( 3 . 0-4 . 0*C0N ( 1 > •I CON ( 2 > ) *D Y 1 

ritind out the values of UpOp T heta and C at every Y-- location but 
t X ” 1.0 i •> e. at the upper edsie o^r the piste 


IF (J.NE.50) GO TO '786 
W R I T E ( 2 5 ;■■ 9 7 ) X >■ T A U p A L P p 1 T R 
F 0 F< M A T ( 2 X p •' X == ' p F 3 . 1 1 X ? '' T 1 M E ~ 
1 V 1 )( p ■' I r R ■■ p 1 2. p / ) 


•■• P F 4 . 2 P 1 X P ■' R E L A X F A C T 0 R 
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WRITE < 25 p 98 ) 




FORMAT ( 6 X p •- K ' p 4 X 

p ' Y 

•' p 5 X ■>> •' U ' P 6 X p ■' y '■ p 7 ) 

<p -'THET 

WRITE( 25 p 93 ) 




F 0 RMAT< 2 Xp 45 < 1 H- 

) ) 



WRITE( 25 p 99 ) (KpY 

( K ) 

p LJ ( K ) p y ( i-( ) p ti-iet ( K : 

5 p CON ( K 

F'ORMAT ( 1 7 p F 6 2 p F 

8 . 4 

■ P F 8 » 4 P F 8 (• 4 ? F 8 . 4 ) 



* 7 / X y ' LI rl ’ ) 


1-0=1 fN) 


DO 


A Li ( K ) =:LI ( !■< 
PU(Kd...-M ) 
CONI IN LIE 
CONTINUE 


) IS AU ( K ) :=U ( K ) f; ATHET < K ) =THET < K ) P AGON ( K ) =CfJN ( K ) 
-U C K ) P F-THET < K p L-f 1 ) -THET ( K) p FtGON < K p L-fl ) =-CON ( K ) 


)lut:Lon has been cslculB’ted out till the UF-F--er lead in si ed«i5e» now 
se instantaneous mean Nusselt and Sherwood numbers are calculated 
ilvinsi the eanr:U (3» 16s) end <3 (■16b) respectively 


NUH~-NUX ( Mi 1 ) 

SHM:=--SHX(Mil) 

DO 132 L==2pMp2 

NUH==NUMi4 . OSKNUX ( L ) i2 . OiiNUX ( Lil ) 

SHH:=SHMi4 . OJlcSHX ( 1... ) i2 <■ 0>!<8HX < Lil ) 

CONTINUE 

B6-l./(3.5rj:iELX) 

NUM--~-NUM:=i<ri6 

SHri:=Si-lMJ|<D6 

kIR I TE < 25 p 1 33 ) NUM p SHM 

FORMAT </p2Xp 'MEAN NUSSELT N0---= ' p F7 .5 .» 2X ? ' MEAN SHERWOOD N0^= 
1 p F7 (. 5 ) 

WR I TE ( 23 ? 889 ) TAU p NUM p SHM 

FORMAT < 3X p F-4 p 2 p 2X p ' ! ' p 7X p F 6 , 4 p 6 K p ■- ! p 6 X p F 6 . 4 ) 

14RITE(25p89) 

F0RMAT(2Xp 5.H !l!-S = ) 


"'hack whe'bher steady state i-s 
1 c u 1 (::• 'L J. o n ‘S * i 'i* r*i o t d o *L o ri e x L 


r e a c h e d o r n o t <■ 
b i. rn e s "L e p- 


i -f 


s t c 


t n e 


DO 990 L=LI. pM 
DO 991 K=:1pN 
A<K)=i.JI (KpL) 

B ( K ) ^=PU ( K p !- > 

CONTINUE 

CALL MAX(ApB.pNpEPS) 

I F (EPS . or . SSMALL ) GO TO 993 
CONTINUE 



. 1 . i; 


c t >[; 1 1 1 1 $ >lc 1 1 -M 1 5{< 'M 5lc t >;■: 1 5}c 5K $ 1C ;{< 1 1 1 ;}: t f. 1 1 1 if; $ ;t: $ :K $ $ :K $ $ if; ;f . :•[; 
S U B R 0 U T I N E T R I D I ( A ? B » C p X p Rl p H ) 


c: 

SOLUTION 

OF' N 

; TIE 

ID It 

JiGONAi... 

TYPE EOUAT 

T r’i i-A a 

.1 iJ i X %.,» 


c: 

A ^ X( 

J--1 ) 

“r Iv« 

$ : 

K ( J ) -f 

c ;{; X(j-fi) 

•■■■’ i\ 

14 HERE 

f' 

K.,' 

A 1 £5 W R J, 

TTEN 

FOIE 

A(. 

,.l ) y h' 

t" OK K K J) ^ 

xf '1* jT-v -v- 1 

J. i„‘ J i 

■!E SOL..!'! VECTEiFo 

c 

ALL THE 

DIMIEN 

!SIO: 

NED 

VARIAB 

LIES l-IAVE D 

I HENS I [ 

■■•A i"!i 0 iE. 1’ •! V 

r'- 

A ( 1 ) SC ( N 

) ARE 

; NO 

I B1 

EFINED 

IN A TRIBI 

riSiONA!,.. 

0../.. 0' 

c 

S R ARE 

hESlR 

;0Y1E 

Dt. 






FUEAL A 

( N ) i? E 

KN) 

p C ( i 

s!) pX<N) 

p R ( N ) p BN 




Ft ( N ) ~mA 

( N ) /E 

! ( N ) 







R ( H ) "-R ( N ) /B ( H ) 
DO 10 I. ""3?N 
Mi 3 ■■•■:[ 

...]•■::!< -I 


BN-i ./(BC J)-'A(K)>i;C( J) ) 

(*■1 ( -.J p "" ri ( -J P iv* N 

1 0 R ( J ) “• C R ( J ) -C ( J ) $R C K) ) i|;BM 

A J . ) K ( 1 ) "" L.- 1 i‘\ \ B \ 1 p H k .' ) ij ( 1 ) ) 

DO 20 :i:==2pN 

2 0 X ( I ) =" R ( I ) - A ( 1 ) X ( I i ) 


FfETURN 

END 


S U B R 0 U T I N IE M A X ( X p Y p H p IE F’ S ) 


DIM ENSIGN 


; ( N ) p ¥ ( N ) 


..J, -J, -I,- 

-T. i*r. /{', 4 ', /f, ,'j'. 


DO 10 I-=1pN 
EPSL^ABS ( X ( I ) ••••¥ ( I ) ) 
EPS==^AMAX1 (EPSpEPSL) 
10 CONTINUE 
RETURN 
END 
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